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We investigate the reconstruction problem of limited angle tomography. Such problems arise naturally in applications 
like digital breast tomosynthesis, dental tomography, electron microscopy etc. Since the acquired tomographic data is 
highly incomplete, the reconstruction problem is severely ill-posed and the traditional reconstruction methods, such 
as filtered backprojection (FBP), do not perform well in such situations. 

To stabilize the reconstruction procedure additional prior knowledge about the unknown object has to be inte- 
grated into the reconstruction process. In this work, we propose the use of the sparse regularization technique in 
combination with curvelets. We argue that this technique gives rise to an edge-preserving reconstruction. Moreover, 
we show that the dimension of the problem can be significantly reduced in the curvelet domain. To this end, we give 
a characterization of the kernel of limited angle Radon transform in terms of curvelets and derive a characterization 
of solutions obtained through curvelet sparse regularization. In numerical experiments, we will present the practical 
relevance of these results. 

Keywords: Radon transform, limited angle tomography, curvelets, sparse regularization, dimensionahty reduction. 



1. Limited angle tomography: Introduction, Organization and Notations 

1.1. Introduction 

Limited angle tomography problems arise naturally in many practical applications, such as digital breast tomosyn- 
thesis, dental tomography, etc. The underlying principle of these imaging techniques consists in two steps: First, the 
data acquisition step, where a few x-ray projections of an object are taken from different view angles (within a limited 
angular range). Second, the reconstruction step, where the attenuation coeflicient / : — > M of the object is approx- 
imately reconstructed from the given projection data. In this paper we are concerned with the second step, i.e., with 
the development of an appropriate (adapted) reconstruction technique which takes into account the special structure 
of the limited angle tomography. 

To this end, we consider the Radon transform as mathematical model for the acquisition process which is defined 



where L(6, s) = |jc e : xi cos 6 + X2sm6 = s^ denotes the line with normal direction (cos 6, sin 6)^ and a signed 
distance from the origin 5 e M. Furthermore, we assume / to lie in the natural domain of the Radon transform, 
i.e., / is such that ([TJ exists for all {6, s). Whenever we write Kg/is) instead of KfiO, s), we consider fifiO, s) as an 
univariate function of the second argument s with a fixed angular parameter 9. In this case, we will call the function 
'Rgf a projection of / at angle 6. 

In contrast to the classical computed tomography, in limited angle tomography the data KfiO, s) is known only 
within a limited angular range, that is, for 6 e [-(t, O] with < n/2. To emphasize that the Radon transform Kf of a 
function is defined only on a limited angle domain [-O, <i>] x M, we will write "Ro/ and call it the limited angle Radon 
transform. As a consequence of the limited angular range, the reconstruction problem y = iRo/ becomes severely 
ill-posed ll22lfTTll . Thus, small measurement errors can cause huge reconstruction errors. 
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This is a serious drawback for practical applications since the acquired data is (to some extent) always corrupted 
by noise. The practical reconstruction problem is therefore given by the equation 

/ = ??<s/ + 7/, (2) 

where denotes the noise, 5 > is the noise level, i.e., ||77|| < 6. The aim is to find an approximation to / from the 
noisy measurements y'\ 

It is well-known that classical reconstruction methods, such as filtered backprojection (FBP), do not perform well 
in such situations meaning that they are sensitive to noise |22|. To stabilize the inversion additional prior knowledge 
about the solution has to be integrated into the reconstruction procedure 1121 . Usually, variational methods are used to 
obtain a regularized solution of the reconstruction problem which is given as a minimizer of the so-called Tikhonov 
type variational functional 

r„(/) = ||^o/-/||2 + «A(/), (3) 

where a > denotes a regularization parameter and A : dom(A) — > [0, oo] is a convex and proper functional ||29l . 
The first term in (|3]l - the data fidelity term - controls the data error, whereas the second term - the so-called penalty 
or prior term - encodes the prior information about the object. 

The choice among the various prior terms and, thus, regularization techniques depends on the specific object 
(which is imaged) and, to some extent, on the desire to preserve or emphasize particular features of the unknown 
object. Usual choices for A are any kind of a smoothness (semi-) norms |29|. For instance, the Besov norm allows 
to adjust the smoothness of the solution at a very fine scale lfT9l l26l l20i . Another prominent example in image 
reconstruction is the total variation (TV) norm which is used in particular for edge-preserving reconstruction, ifTJlfTTl . 

Indeed, to preserve edges is an important issue for medical imaging. However, it was pointed out in fTSl . that TV 
reconstruction may be not an appropriate choice for medical imaging purposes. One reason for this is that TV regu- 
larization favors piecewise constant functions and, hence, produces staircase effects (cf. ['27', "81) which may destroy 
relevant information. Hence, piecewise constant functions may be not appropriate for our purpose. To overcome this 
problem, higher order total variation priors were considered by some authors, see for example 1 1 1. In this work will 
use curvelets to avoid such problems while preserving edges of the reconstruction. 

Another issue we are concerned with is the fact that in limited angle tomography one can not expect to get a perfect 
reconstruction (though the limited angle problem is uniquely solvable in some mathematical settings). Depending on 
the available angular range some structures of the unknown object can be reconstructed (are visible) and some can 
not be reconstructed (are invisible) ||23]| . To our knowledge the information about the visible and invisible structures 
(which is encoded in the data set) is not exploited by any of the mentioned reconstruction methods. 

In view of the above discussion, our goal in this work is to design a reconstruction method for limited angle 
tomography which is 

(i) stable, i.e., insensitive to noise, 

(ii) independent of acquisition geometry, 

(iii) edge-preserving, 

(iv) adapted to the limited angle setting, i.e., exploits information about visible and invisible structures. 

First thoughts on this topic have been formulated in an extended abstract (2 pages) that is submitted to the Proceedings 
in Applied Mathematics and Mechanics, U5J . 

1.2. Organization of this paper 

In the first part of Section|2]a brief description of the curvelet dictionary will be given. The second part of Section 
|2]the technique of sparse regularization will be introduced as a stable reconstruction method. In Section [3] we will 
discuss the relation between the curvelet sparse regularization technique and the curvelet thresholding which was 
proposed in |4 |. Based on this discussion, a characterization of the curvelet sparse regularization will be given for the 
full angular problem by using the biorthogonal curvelet decomposition (BCD) [4i. Afterwards, the limitations of the 
BCD approach will be discussed. 
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Our main results will be presented in Section |4] Here, we will first prove a characterization of the kernel of 
the limited angle Radon transform in terms of curvelets. As a consequence, a characterization of curvelet sparse 
regularizations will be derived. These results will be applied to a finite dimensional reconstruction problem in Section 
|5] By performing dimensionality reduction of the reconstruction problem in the curvelet domain, an adapted curvelet 
sparse regularization approach will be introduced. In Section|6]we will discuss some of our results. 

We will conclude this paper by showing some numerical experiments in Section |7] In particular, we will show 
that the execution times of the adapted curvelet sparse regularization significantly reduces while preserving the recon- 
struction quality. 

1.3. Notation 

We state here some notations which will be used throughout this paper: 

The inner product of x,y e W will be denoted as jc ■ y or simply xy. When not otherwise stated, inner product in a 
function space X will be denoted by (/, g}x. The norm of a vector jc e M" will be denoted by |x| whereas the norm in 
a function space X will be denoted by 

We will be using some classical function spaces, such as the space of Schwartz functions >S(M") and the spaces of 
measurable functions LP{D), without reference since they can be found in every book on functional analysis. Same 
holds for the classical sequence spaces 1^. 

The Fourier transform / of a function / e iS(M") is defined by 

m^(2n)-"^^ r /«e-''-fdx. 

The inverse Fourier transform is given by fix) - fi-x). Basic properties of the Fourier transform will be used without 
proof. For details about the Fourier transform we refer to ll30ll . 

For 77 e [0, 2ti\, we define to be the rotation operator g,jf{x) = f{R,jx), where the rotation matrix R,j is defined 

by 

^ _ / cos T] sin7/\ 
'' \^-sin/7 cosT/j' 

Eventually, we refer to f22] for notations and some basic facts about the Radon transform. 



2. Stabilization of limited angle reconstructions by sparsity in the curvelet domain 



In this section we are going to address our goals (i)|(iii) stated at the end of Subsection |1.1| To stabilize the 



inversion we need to incorporate some a priori information into the reconstruction which permits an edge-preserving 
reconstruction, or at least, does not smoothes edges in the reconstruction. We will do so by assuming that the functions 
we are going to reconstruct belong to a class fi^ which consists of functions that are except from discontinuities 
along curves. To translate this qualitative information into a mathematical language we use the fact that functions 
in £^ are optimally sparse with respect to the curvelet frame Q. Hence, the technique of sparse regularization ifTOl 
seems to be appropriate in our setting. 

To this end, we briefly recall the definition the curvelet frame [7J and collect some basic facts about technique of 
sparse regularization. 

2.1. The curvelet dictionary 

At scale 2"^, j e No, we first define the generating curvelets i^j.o.o^ in the frequency domain using polar coordinates 
(r, (jj) by 

_ /2r;72i+i \ 

.A,;o,o(r, u>) = 2-3-'/4 . W{2-J ■ r) ■ y w , (4) 



where W{r) is a radial window and y(a») denotes an angular window. The windows W and V are both real and smooth, 
i.e., W, y e C°°. Furthermore, we assume that supp W c (1 /2, 2), supp V c (-1, 1) and that the following admissibility 
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conditions are satisfied, 



W^iVr) = 1, re (3/4,3/2); 

oo 

^y\w-0 = l, we (-1/2, 1/2). 



/=-oo 




Figure 1: Support of curvelets in the Fourier domain for j = I (dark gray), j = 3 (gray) and / = 5 (light gray). 
The family of curvelets {lAj,/,*:} now constructed by translation and rotation of generating curvelets i/'/,o,o- That 

V J jj,k 

is, at scale 2"^, the curvelet i/ry / j; is defined via 



(5) 



where Rg.^ denotes the rotation matrix (cf. Section 1.3 1 with respect to the scale-dependent rotation angles Ojj and 
scale-dependent locations which are define by 



l-n-2 



<ijm+i 



<l<2 



r>/2i+i 



k = {h,k2) e I? 



Since the window functions W and V are compactly supported, and in particular, since the support of W{2^-) is 
contained in (1/2, oo), it follows from (j4| and (|5]l that, in the Fourier domain, each curvelet is supported on a polar 
wedge which has a positive distance to the origin, see Figure [T] We have ^j^kiO = for all |^| < 1 /2 and for all 
admissible indices (j, I, k), i.e., the region Ucy.;,*) supp ijjjiji covers not all of the M^. Thus, the system {lA ;,/,*■} does not 
contain any low-pass element. 

To complete the definition of the curvelet system we define the generating low-pass function i/'-i,o,o in the Fourier 
domain by 

oo 

iA-i,o.o(r, u) = Wo(r), wl{r) := 1 - ^'(2"-''-) 
and complete the curvelet system by all of its translates o./tj^-g^^- 



Remark. 

In the spatial domain, the essential support of curvelets is an ellipse which is located near bf and oriented along 
the orthogonal direction 0^^ - Ojj + n/l. The directional localization becomes higher when the scale parameter j 
increases. Thus, curvelets are highly oriented at fine scales, see Figure|2] 
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(a) lAsoo 



(b) i/'fiso 



Figure 2: Curvelets at different scales and different orientations. Left image shows a curvelet with orientation 6j o = 0° whereas 
the right image shows a curvelet with orientation 06.5 = 56.25°. 

The index set of the completed curvelet system is now given by 

l^[(-l,0,k): keZ^]u{(j,l,k): j eNo, k eZ', -2^^'^^^^ < I <2^j'^^^^] -. lo U Jj. (6) 

Note that each index (j, l,k) el has a 3 parameter structure, where j denotes the sca/e-parameter, k = (^1,^2) is the 
location parameter and I is the orientation parameter The system {(Aj,/,*}!, , j^^j is now complete in the sense that it 
constitutes a tight frame for L^{M?), |7 1. For each / e L^(]S?) there is a curvelet representation 

/= Z {>Pjj,kj)>1fjXk (7) 

and a Parseval relation holds, 

U,l,k)el 

We conclude this section by noting that curvelets can be understood as further development of wavelets ||9l . 
Therefore, as it is well known from the theory of wavelets, the curvelet dictionary provides a sparse representation of a 
large class of functions. This property qualifies curvelets for the use within the framework of the sparse regularization. 
An even more important property of curvelets lies in the fact that they offer an optimally sparse representation of 
functions that are except form discontinuities along curves (5). This means, that curvelets encode edges in a 
very efficient way. Thus, sparse representations of functions with respect to the curvelet dictionary can be considered 
to be edge-preserving. 

2.2. Curvelet sparse regularization ( CSR) 

Our aim is to solve the limited angle reconstruction problem (|2| in a stable way such that the edges are preserved. 
To this end, we use sparse regularization of curvelet coefficients. The idea of sparse regularization is to determine a 
solution / of the problem which is sparse or compressible with respect to the curvelet frame. Sparsity of / means that 
the series expansion ([7]l of / has only a very small number of curvelet coefficients which are non-zero, whereas 
compressibility of / means that / can be well approximated using a sparse series expansion. 

To this end, we have to formulate the reconstruction problem in the curvelet domain, i.e., we are interested in the 
recovery of the curvelet coefficients cjj^k = {^jj,k,f) of / instead of the function / itself. We assume that the unknown 
object / can be represented by a finite linear combinations of curvelets, i.e., / = Tin=i i'Pn,f) lAn with n - n{j, I, k). 
Further, we define the analysis operator T and the synthesis operator T* as 

N 
«=1 
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Then, the reconstruction problem (|2]) can expressed in terms of curvelet coefficients via 

y^^Kc + rj, K-.^'R^T*. (8) 

A solution to ([HJ by sparse regularization of curvelet coefficients is given as a minimizer of the -penalized 
Tikhonov type functional, i.e., 

,1,2 



c = arg min - \\Kc - y\,^,,^^^ + ||c||i,, , (9) 

where ||c||i,^, - 'Yjk'^k\ck\ denotes the weighted 1-norm with a weight sequence w satisfying Wk > wo > 0. A 
reconstruction for the original problem Q is ffien given by applying the synthesis operator to the regularized curvelet 
coefficients c, i.e., 

N 

f^Yj^c„ilf„. (10) 

n=I 

We note that sparse regularization is indeed a regularization method lfT0l . ||29l Sec. 3.3]. Therefore, the computation 
of a reconstruction by (|9]) and ([T0| is stable and favors sparse solutions lfT0l[T6ll . We will refer to this method by the 
term curvelet sparse regularization or CSR, respectively. In the previous subsection, we have discussed that sparse 
representation of functions with respect to the curvelet dictionary are edge-preserving. Consequently, curvelet sparse 
regularization gives rise to an edge-preserving reconstruction method. 

In the following proposition we give a general characterization of minimizers of the /"'-penalized Tikhonov func- 
tional. Though the proof can be found in |16|, we will recall it here for the sake of completeness. To this end, we 
define the so-called soft-thresholding operator ■ — > by 

(S„(x))k = S^i^ixk) := max{0, M - Wk] sgn(xk). (11) 

Proposition 2.1. The set of minimizers of the (^-penalized Tikhonov functional 

is non-empty. Furthermore, each minimizer cof^ is characterized by 

c^Sy,.{c-yK*(Kc-/)) (12) 

for any y > 0. 

Proof. We follow the proof of lfT6l . Since *P is convex and coercive it follows that there is a minimizer c of *P. We 
denote by df(c) the subdifferential of a convex function / at x. Each minimizer c is characterized by the requirement 

(cf. EH]) 

QedW)-K*iKc-g) + d\\c\h,, 

which is equivalent to 

-K*(Kc-g)ed\\c\U,,. 

Multiplying by y > and adding c to both sides yields 

c - yK*(Kc -g)ec + r-5||c||i,, = (id + yd\\-\U^„)c. 

Following the arguments in (16] we get that (id + yd |Mli,vy)"' exists and is single valued. Hence, the above inclusion 
is characterized by the equation 

c = (id + yd \\-\UJ-\c - yK*(Kc - g)). 



A simple calculation shows that (id -H yd |Mli_„.) ' - S 



yw 
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3. Relation to biorthogonal curvelet decomposition (BCD) for the Radon transform 



In this section we will show that if the data is available from the full angular range, i.e., if we are dealing with 
the Radon transform H rather than the limited angle Radon transform H,^, an explicit formula for the minimizer of 
the /"'-penalized Tikhonov functional (|9]) can be derived using the biorthogonal curvelet decomposition (BCD) for the 
Radon transform |4|. This formula is closely related to the BCD based reconstruction which was also proposed in Q. 
Afterwards we will discuss that the curvelet sparse regularization can be understood as a natural generalization of the 
BCD reconstruction. 

3.1. Full angular range 

We now briefly recall the definition of the BCD for the Radon transform. For details we refer to pl. If not 
otherwise stated, we let n - n(j, l,k) el and denote the curvelet frame by {ifr„}. In order to derive the BCD for the 
Radon transform a pair of frames {U„} and {V„) is constructed for ran CR) c L-(Mx5') such that 

and a quasi-biorthogonal relation {V,,, ^/;^')L-(Rxs') = 2^"^' {^n,^n'), holds for all «,«' e J. In particular, there is an 
L^-norm equivalence property 

neJ 

for all g e ran {H). Similar relations hold for {V„]. Using these notations, the BCD of the Radon transform is given by 
the following reproducing formula 

/ = 22'<^/'^«>iW)'An, (13) 
neJ 

where / € L^(1R^) is assumed to be a finite sum of curvelets {i/f,,} lH. Note that the curvelet coefficients of / are 
computed from the Radon transform data Hf. 

We now use the above the frames (?/„), (V,,) and its relations with H and respectively, to compute the minimizer 
of the /'-penalized Tikhonov functional. We assume that e ran('??) and denote (■,■) - (■,-)^2(k2) and [■,■] - 
(; ■)l^(Rxs|)- Withy^ = [/, U„l we get 



nel 

= Yj\^'Rf,U„]-\y',Unf 

nel 

= ^|</,!R*f/„)-[/,f/„]|' 

= ^|(/,2-Vn)-[y',f/«]| 

nel 



nel 



Using the definition of ||-||i_„, we see that 

nel 

with a suitably chosen constant a > 0. Hence, we have 

c = ai-g min ||/:c -/||^^ + \\c\\^^„ 



arg min V (|2 ■'c,, - y^f + aw„ |c„|) , (14) 

ceW I 
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which can be minimized by minimizing each term in ( [T4| i separately. Note that each term in ( [T4] i is of the form 
\ax — bf + c\x\ and its minimum is given by S c/(2(P-}(b/a) where S c/(2a-) is the soft-thresholding function from ( |TT] l 
with the threshold c/(2a^). Therefore, the sparsity regularized curvelet coefficients c are given by 

We have now proven the following Theorem. 

Theorem 3.1. The solution of the full angular problem — "Rf + rj via curvelet sparse regularization is given by the 
( closed) formula 

f^Y,^2V-.,,JVyl)^n. (15) 

nel 

The relation between the reproducing formula ( |T3] l and ( [T5] l is now obvious. If the thresholding parameters 
'Wn - w„{6) in ([TSj are chosen such that w„{S) — > as 5 ^ 0, i.e., they vanish if there is no noise present in the data, 
then ( [T5| ) reduces to ( [T3] l. On the other hand, if the data is corrupted by noise, then, the curvelet sparse regularized 
solution is simply a thresholded version of the BCD reproducing formula. The stabilizing character of the curvelet 
sparse regularization is reflected by the inherent thresholding of the curvelet coefficients (see also ([T2ji). 

In |4J a very similar reconstruction rule was derived. Starting form the BCD reproducing formula the authors 
proposed to use soft-thresholding of coefficients in ([T3|, i.e., 

nel 

with a scale dependent threshold tj. We see that this formula coincides with ( [T5] l for a suitably chosen thresholding 
sequence r = (jj). 

Remark. 

Note the ill-posed nature of the reproducing formula ( [T3| l. This is evident because the coefficients 2-' {llf, f^/j)^2(Rxsi) 
corresponding to fine scales (large f) are amplified by the factor 2'. Since noise is a fine scale phenomenon, there will 
be very large reconstruction errors when the data is corrupted by noise. 

3.2. Limited angular range &- Limitations of the BCD 

We have seen that there is an explicit expression ( fTS) for the CSR reconstruction in the case of full angular 
tomography. We also noted that the thresholded BCD reconstruction ([T6| leads (under certain conditions) to the same 
reconstruction. To extend this observation to the limited angle tomography a biorthogonal curvelet decomposition for 
the limited angle Radon transform would be needed. To our knowledge there is no such BCD available for the limited 
angle Radon transform. Consequently, in the case of limited angle tomography, the CSR reconstruction ( [TO] i can not 
be expressed explicitly as it was done for the full angular range in ( [T5] l and the BCD reconstruction of Candes and 
Donohod 14] can not be applied in this situation. 

In contrast to the BCD method, a reconstruction of the limited angle problem can be computed using CSR. Hence, 
curvelet sparse regularization can be understood as the natural generalization of the thresholded BCD reconstruction. 

Curvelet sparse regularization offers even more flexibility compared to the BCD method. For example, the imple- 
mentation of the thresholded BCD method is difficult for acquisition geometries which are different from the parallel 
geometry. This is because the BCD method requires discretization of the functions U„ which live in the Radon do- 
main. The implementation of the curvelet sparse regularization approach, however, is independent of the acquisition 
geometry. One needs only to implement the system matrix. Moreover, the generalization to higher dimensions is also 
easier accessible via curvelet sparse regularization approach. 



4. Characterization of Umited angle Radon transform 

In Section|2]we presented curvelet sparse regularization as our method of choice for the limited angle tomography 
because it is stable, edge-preserving and flexible. In Proposition 2.1 we noted the existence of a solution and showed 
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that each minimizer of the /"'-penalized Tikhonov functional (|9| is given as a fixed point of some operator (cf. ([T2]l). 
This characterization is generic in the sense that it does not take into account the special structure of the underlying 
problem. 

The goal of this section is to give a characterization of the minimizer (|9]l which is adapted to the setting of limited 
angle geometry. In the following we will show that, depending on the available angular range, a big portion of the 
curvelet coefRients of the CSR reconstruction are zero. 

We state our main results first and postpone the proofs to the end of this section. 

Theorem 4.1. Let < <1> < 7r/2. We define the polar wedge Wtj, by 

Wo = e : ^ = r(cos w, sin co), r e M, |w| < <d} . (17) 
Moreover, we define a proper subset of the curvelet index set by 

^-invisible ^ |(^. IJ^-^^J, supp 0- n Wo = 0} , (18) 



where ipjXk denotes a curvelet and I is the curvelet index set (cf. Subsection 2.1 ). Then, 

5Ro<Ai,a = Ofor all (J, I, k) e J■™^'"^ (19) 

The above theorem characterizes a subspace of the kernel of the limited angle Radon transform in terms of 
curve lets. Using this Theorem 4.1 a characterization of curvelet sparse regularized solutions to the limited angle 
problem can be derived. 

Theorem 4.2. Lef < <t< n/2, y^ e ran (-Ro) and let jjjJ^^iW': be defined by Then, 

1 .||2 



c = argmin ILzf^.^R, + IWI,,., 

satisfies 

cjj,k^Oforall (j,^,Oe J|^;™™^ 



We start to develop the proof of Theorem 4. 1 first. To this end, we need some auxiliary results. Though the content 
of the following lemma is classical we will give a proof for the sake of completeness. 

Lemma 4.3. Let b > Q and (5* be a function defined by 



27r Xi^ 



5''(x)= - ( e-"-^d^. (20) 



Then, it holds that 5 — > 5 pointwise in »S'(M) as b ^ oo, i.e., for ip e S(M.) we have 



Proof. A simple calculation shows that 



lim 6''(x)ifi(x)dx = ifiiO). (21) 



Let (fi e S(R). We split the integral in ( |2T] i 

1 psm(/«) 1 rsm(to) If sinibx) 

— I ip{x)ax— - (p(x)ax-\ — | ip(x)ax 



=: I lib) + hib). 
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Observe that g(x) - X{\-\\>£}ix)f(^)/^ is in i'(II^) and it follows from the Riemann-Lebesgue lemma ll30l Theorem 
1.1 + 1.2] that 

1 r [2 
lim hib) - lim - I sm(bx)g(x)dx - \ - lim Im(f(fe)) = 0. 

To compute limi^co 1\{V) we use the Taylor expansion of at x = 0, ip(x) - i^(0) + '(^^x with x e (-e, s) and 
e [-.X, x]. Note that ip' is bounded and, hence, (p'i^^) is integrable on (-e, e). 



lim hib) - lim 




x)^(x) dx 


= lim j 




x)(^(x) + (;o(-x)) 









_ 2^(0) 


lim 


r sin(/7x) 

( dx + 


TT 




io 


2^(0) 


lim 


r"-^ sin(x) 


n 




Jo 









where we have again applied the Riemann-Lebesgue lemma to the function ^(x) = X(%e){x){ip'{^x) + f'{-^x)) and have 
used the asymptotics for the sine integral, limj^oo sin f /f df -njl. ■ 



The key observation for the proof of Theorem |4.1| is contained in the following lemma. 
Lemma 4.4. For f e kS(]R^) and £, - (cos 77, sin 77), - (- sin cos r]), we have 

£f(t-^)dt^ Jj(t-t)dt. (22) 

That is, integration in the spatial domain along a line through the origin corresponds to the integration along a 
perpendicular line through the origin in the frequency domain. 



Proof. We first show that 



r /(xi,0)dxi = r /(0,^2)d^2. (23) 

Jk Jr 



Using the notation S''{x) - (In) ' j^^ e'"^ d^ we compute 



X 



1 



/(Xl,0)dx, = ^ I I /(^)e"'^'d^dX; 



lim r r /(^)e"'f> d^dx, 

In b^ca J_Ij Jjj2 



= lim /(^i,6)5''(^i)d^idf2 



lim /(fi,ft)5''(^i)d^, d^2 



/(0,^2)d^2 



where we have used Lemma 4.3 In the above computation, the change of the integration order and the interchange 
of the limit process and integration are allowed due to Fubini's theorem and the dominated convergence theorem, 
respectively. 
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Since Fourier transform and rotation commute, we get using notations in S ubsection 1 1 . 3 1 and ( |23] l that 

X/"«'"=//h-(o))'" 

Jr 
Jr 

Jr 



Using Lemma 4.4 we are now able to derive a formula for the Radon transform of curvelets. 
Theorem 4.5. Let ij/jj^k be a curvelet (cf. ^) and denote 0(a>) — (cos co, sin w)'''. Then, 

/orj/zi+i \ _ 

Kifrjj^tiOico), s) = — - — (w + 9jj)\ yl2^w[V {bf, 6{oj + 9jj)) - 2's) , (24) 



where b'j^' and 9j i are defined in Section 



2.1 



Proof. Let 6 6(cj) - (cos o), sin coy, a> e [-<I>, O], and Tpf :=/(• + p), p e M^. First note that each curvelet t//jjii is 
a Schwartz function since, per definition, its Fourier transform is C°° and compactly supported. Hence, we may apply 
Lemma 1441 

'R^i/fjj,k(0,s)^ f <Ay,a(*0 + f0^)df = f(T,,eif^jj,k)(t0^)dt^ f e'<''''"'>iAy.«(f0)df. (25) 
Jr Jr Jr 

At scale 2"-', each curvelet ifrp^k is defined via translation and rotation of a generating curvelet i/'/o^o, cf (|5]l. Using 
the relation of the Fourier transform and rotation as well as the relation Tpf(^) = e''^P'^\f(^) we see that 

^^jjA^) = e"'<'-*''^'^V,-,o,o(«fl,/). (26) 
Now plugging ( p6l l into (|25|, together with Q, we deduce 

'R^^j.i,ki0,s) = r e''V-'(''-%'<"")iA;,o.o(^fl,,(f0))df, 
Jr 

= r e''"e-''(''^''«("+'''''>)./r;,o,()(r, + %/) df, 
Jr 

= 2-3-''^V I ^^^(w + 0,;,) I J •»("+«^')>lV(2--'f) df. 



^2r>/2i+i 

^—iaj + 0jj)\ J e-''P'(^-''("+»^.')>-V(r)dr, 

V2^ (2^' {bf, 0icj + 0jj)) - 2h) . 



The proof Theorem 4. 1 is now a simple consequence of Theorem 4.5 
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Proof of Theorem l4.1l Let j e N. To abbreviate the notation we let aj - n^^l^j^^^ and denote by 

[-7T, -;7r + O] U [-O, O] U [tt - O, n] 

the symmetric (visible) angular range of the limited angle Radon transform "Ro (cf. Figure|3]l. 

According to Theorem 4.5 we have to determine all Op such that V{aj(- + 0j,/))Uo = 0. Since suppV c (-1, 1) 
and Ojj e [-tt, tt], we have 

0jj i (-aji - aj' - w) => V(aj( (jj + 0jj)) = 0. 

for all (L) e Ao- Therefore, by defining (cf. Figure [3]l 

Aoj := [-7r,-7r + (O + flj')] U [-(O + aj'), O + aj'] U [;r - (O + flj'),7r], 

we see that V{aj(- + 0j;/))U« = holds whenever 6jj i Aoj. The assertion follows by defining the invisible index set 
of curvelet coeflicients as 



-TT — TT + tE> 

-i 1 



—I— 



- * 



Figure 3: The symmetric (visible) angular range of the limited angle Radon transform, A^, and its scale-dependent version, A^ j. 



We summarize the steps that were needed to prove Theorem |4. 1| This procedure is also illustrated in FigurefflTo 
evaluate "Roi/';,',it(^' for a curvelet lA;,;,^, the integration was shifted to the Fourier domain according to Lemma |4~4| 
This related the value 'Rm^ jjjS.^, s) to the integration of i^^ / ^ along the line - {f^^ : t e M). Because of the limited 
angular range, the union of all such lines, 

(27) 



u 



j;e[-'I),<I)] 

covers not all of the M^. Above, we have again used the notation ^{rj) - (cos t;, sin 77). To prove the assertion, we 
computed all those curvelet indices ( j, I, k) such that supp t/fjj^t n Wo - 0. 



Now, we turn the proof of Theorem 4.2 This will be a simple consequence of Theorem 4.1 and the following 
lemma. 



Lemma 4.6. Let f : [ —00, 00] be defined by f(x) — 2«eJ fi^n), where ip 

that f is proper. Then, it holds that 

y e df(x) o y„ e dip(x„)fijr all n e I. 



. is a convex function such 



(28) 



Proof. First note that, since / is proper, for x e {^{I) we have df(x) = if f(x) = 00. In what follows we therefore 
assume without loss of generality that f(x) < 00. 

Suppose y„ e dtp(x„) for all n e I. Then, by definition of the subgradient we have 

Vz„ € M : (p(z„) > (fiixj + y„(z„ - x„). 

Summing over n implies 

Vz e : ^ ifi(zn) > ^ (ip{x„) + y„(z„ - x„)) = ^ (p(x„) + {y,z-x). 



nel 



nel 



nel 
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which is by definition of / equivalent to y e df(x). This proves the implication of the statement. 

On the other hand, if y e df{x), then f(z) > fix) + {y,z- x) for all z e (^(I)- In particular, this holds for aU 
z = X + he„ with /z e M and n e I, where e„ - (6i,„)iei and 5, „ denotes the Kronecker delta. Therefore we have 

V/i € M Vn € J : fix + he„) > fix) + {y, he„) 

o ^ ipixj + 6i,„h) - ^ ipixi) > y„h 

iel iel 
o ipix„ + li) - (pix„) > y„h 
<=> yn e d(pix„). 



Proof of Theorem 14.21 As in the proof of Proposition 2. 1 we see that c fulfills the following relation 



-K*iKc-/)ed\\c\U,,.. (29) 
Since y^ e ran ('??)j, we have that y^ - Kc^ = 'R^T*c^ for some curvelet coeflicient vector e M'^. Thus, 



-K'iKc - /) = -K'Kic - c^) = -K' 



^ic- c'^)jik'RQ>tf/j^, 



By Theorem 4. 1 it holds that 



xjj^k = for all (;, k, I) e JJ"-^""". 



(30) 



Now let ij,l,k) e Ji^v'^*ib^^ in view of ((29| and ((30]l, Lemma 4.6 impHes that e d[wjii, \cjj^k\)- However, this 
means that cjj^k minimizes the function wp^k \ ■ \ and since Wjj,k > we get that cjj^k - 0. m 



5. Adapted curvelet sparse regularization 

In this section we are going to apply the results from Section[4]to a finite dimensional reconstruction problem. We 
will show that, in this setting, a significant dimensionality reduction can be performed. Based on this approach, we 
wiU formulate the adapted curvelet sparse regularization (A-CSR). 
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5.1. Discrete reconstruction problem 

We consider the discrete version of the reconstruction problem. To this end, we model / as a finite linear combina- 
tion of curvelets, i.e., / = Yjn=i '^n^n, where n = n(j, I, k) is an enumeration of the curvelet index set I (cf. Subsection 



2.1 1. Moreover, we assume to be given a finite number of measurements y,„ = 'Rt^iOm, Sm), 1 < m < M G N. Then, 



each measurement y,„ can be expressed as 

N 

ym = ^o/(fm, S„) = ^ CnK^il/nifim, ^m)- (31) 
(1=1 

Now, let us define the so-called system matrix K by K„, n - 'Rq,il/„(6m, s,„) for 1 < m < M and n e I. Then, the discrete 
version of the limited angle problem (|8]l reads 

y^Kc + rj. (32) 

Note the abuse of notation. In contrast to ([8|, where K - HcsfT* denotes a continuous operator, here, K is its discrete 
version. We want to point out, that the reconstruction problem (|32]) is formulated in terms of all curvelet coefficients 
c = T f. That is, to solve ( (32] l, we need to compute c„(jj^k) for all possible curvelet indices {], I, k) € I. The dimension 
of the reconstruction problem, given by = \I\, does not depend on the available angular range. In what follows, we 
will use the method curvelet sparse regularization to solve this problem. 

5.2. Dimensionality reduction Cr Adapted curvelet sparse regularization (A-CSR) 

First, note that the results from Section]?] are formulated only in terms of the angular range parameter In turn, 
this parameter is completely determined by the acquisition geometry. Hence, it is known prior to the reconstruction 
and can be extracted from the given data by 

O - min {(fi : 3(po € [-n, tt] : V 1 < m < M : 0„, e [(po - ip,(pQ + cp]] . 



Knowing O, we can use Theorem 4.2 to identify those curvelets which lie in the kernel of the limited angle Radon 



transform "Rc^. Their index set can be precomputed according to {I8f or, equivalently, by 

jinvisible ^ |(^. 11^^^ J, (cos 0jj, smOjjf i IVoj} , 

where W^j' - e : ^ = r(cos w, sin iiSf , r e M, |w| < <1> -H nT ^■'^^^ } is a polar wedge at scale 1^' and Qjj is 
the orientation of the curvelet ^ j^ik- In what follows, curvelets i/^j/^ as well as curvelet coefficients c„^jjj.-, with 
( j, I, k) € J^jj^i'^'bie ^jjj called /nv/i/fe/^from the given angular range. Accordingly, the index set of visible curvelet 
coefficients is defined by 

7-visible 7- \ 7-invisible 



In view of Theorem 4.1 it holds that K,„^„ = 'R^iff„(e„„ s,„) = for n e jmvisibie ^j^^ f^j. all 1 < m < M, i.e., 
those columns of the system matrix K which correspond to the invisible index set are identified to be actually zero. 
Therefore, we may define a new system matrix /To with respect to the visible index set by 

(K^)„r,„ = n^Mdm, sj, 1 < m < M, « e J™"''^ 

Such a reduced system matrix Kq, has the size M x Since |JJ[I^*''^| - \I\ - \r^^^*'^'^\, the number of columns 

is reduced by the quantity Using the reduced system matrix we formulate the adapted (or reduced) limited 

angle problem as 

/ = K^c + 7]. (33) 

The dimension of the adapted problem, given by A^o = l-^l - li^JU^'^*'*^!, now depends on the angular range parameter 
O. From the definition of J^i'^'bie it is clear that as the angular range becomes larger the number of visible curvelets 
increases. Hence, the dimension of the adapted problem N<:, increases as the angular range increases and vice versa. 



'We adapted the term invisible from [241. 
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Apply the technique of curvelet sparse regularization (cf. Subsection 2.2 1 to the reduced problem ( [33| l we formu- 
late the adapted curvelet sparse regularization (A-CSR) as 



Co = arg mm 



\\\K^c- 



l|c|li,„. 



(A-CSR) 



Apparently, the computational amount decreases by using the A-CSR framework instead of the CSR. However, 
according to Theorem 4.2 the reconstruction quality is preserved. In Section]?] we will present some practical experi- 
ments concerning these issues. 



Remark. 

Note that the characterization of Theorem 
index of summation I by J™'''''^. This yields a closed form solution for ( |A-CSR| l 



4.2 



may also be applied to the closed form formula ([TSj by replacing the 



6. Discussion 

This section is devoted to the discussion of the results which were presented in the previous section as well as 
their implications for the practical application of the curvelet sparse regularization. 

General angular ranges. So far, we have worked with a symmetric angular range [-O, <1>] with < C) < njl which 
was centered at <I)o - 0. The results of Section]?] however, can be easily adapted to a more general situation, 
where the available angular range [Oq - O, <t)o + "t] is centered around an angle Oo e [-n, n] . To this end, let r<p„ 
be the translation operator defined by T(^^'R(^ f{9, s) - !Ro/(S + 4*0, s)- Then, the limited angle Radon transform 



with respect to a general angular range [Oq - O, -h <I>] x M is given by Tafgli^. Theorem 4.1 and Theorem 4.2 
can be now applied to T^^'Ri^,, yielding a general index set of invisible curvelet coefficients 

jinvisible ^ 1^^. ; e J . ^^pp ^ ^ /?<Po Wo = ©} , 

where /?o„Wo - {R<s>o^ '■ ^ e Wo) is a rotated version of Wo. 



Computation of the system matrix. In Theorem 4.5 we have derived an expression for the Radon transform of a 
curvelet This expression can be used to compute the entries of the system matrix O analytically, if both, 
the angular window V and the Fourier transform of the radial window W are known analytically. This is useful 
for practical application since, in this case, the system matrix can be precomputed and needs not to be set up 
in every iteration of the minimization of the ^'-penalized Tikhonov functional. This may yield an additional 
speedup of the algorithm. 

Additional stabilization of the limited angle problem. Adapting the problem to the limited angular range has an 
additionally stabilizing effect. This comes from the fact that the reconstruction problem (]33]l is formulated with 
respect to visible curvelet coefficients only. In this way, a big portion of the null space of the system matrix 
(limited angle Radon transform) is excluded from the formulation of the limited angle problem. Therefore, the 
condition number of the reduced system matrix improves which induces an additional stability. 

We want to point out that the formulation of the adapted problem (J33]l does only depend on null space analysis 



of the limited angle Radon transform in terms of curvelets (cf. Theorem 4.1 1. Thus, the adapted limited angle 
problem (]33]l is not related to any reconstruction algorithm. Therefore, the additional stabilization will be 
present if any other method would be used for solving (J33|. The adapted formulation of the reconstruction 
problem ((33]l can be therefore interpreted as preconditioning procedure. 



Related work. In lfl4ll . the adapted curvelet sparse regularization was introduced by using microlocal analysis. There, 
a qualitative characterization of visible curvelets was derived from the characterization of visible singularities 
of E. T. Quinto |23 1 and the "the resolution of wavefront set property" of the continuous curvelet transform ||6l. 
These results were stated there without proofs. 
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7. Numerical experiments 



This section is devoted to the illustration of our results which were presented in the previous sections. To this end, 
two types of numerical experiments were made. In the first part of our experiments we will illustrate the visibility of 
curvelets under the limited angle Radon transform and show how this leads to a dimensionality reduction in the limited 
angle reconstruction problem. In particular, these experiments are meant to illustrate Theorem |4.1| and Theorem |4.2| 
The second part of our experiments is devoted to reconstructions via CSR, A-CSR and filtered backprojection (FBP). 
A comparison of these reconstructions will be presented in terms of execution times and reconstruction quality. 

7.1. Implementation of the minimization algorithm 

For the minimization of the /''-penalized Tikhonov functional (|9]) we implemented a variant of the well known 
iterative soft-thresholding algorithm ifTOl l2l. This algorithm is given as a fixed point iteration of the equation ( [T2| i, 
namely 

c«+i = Sr.. [c" - s„K*(Kc" - /)) , (34) 

where we have used c" s as an initial guess. This procedure consists of a gradient descent step with a subsequent 
soft-thresholding with respect to the sequence t" = (j"ji]^^)(j,i,k)ei- The step length s„ of the gradient descent step was 
chosen such that Q<s_<s„<l<2l \\K\^, [2J. Usually, the thresholding sequence r" is chosen as t" - s„Q w, 
where O denotes pointwise multiplication of the step length i„ and the /"'-norm weight sequence w (cf. Q). This 
weight sequence is a free parameter and has to be selected appropriately because it affects the reconstruction quality. 
In general, there is no rule how to select such a weight sequence. In practice, this often done by trial and error. 

In our implementation we got rid of this weight sequence by choosing the thresholding sequence r„ adaptively 
and scale-dependent at each iteration n via 

=2^'^--^>/V ^2 log, A^„, (35) 

where cr denotes the standard deviation of the noise rj, Np is the number of curvelet coefficients at scale 2"^ and at 
orientation Ojj and 7 e N is the largest available scale parameter for the image size of interest. This thresholding 
strategy was initially presented in |4, Sec. 6]. Since it is based on the white noise model, we assumed throughout our 
experiments the noise to be white Gaussian. 

Moreover, we simulated a practical situation by assuming that the noise, and especially its standard deviation cr, 
is not known. In order to automatize the reconstruction procedure, we used the median absolute value (MAD) to 
estimate cr (cf. ||2T1 p. 565]) by 

cr~ 1.4826 ■MAD(c'j). 

Above, MAD(c'j) is the median of the absolute values of the curvelet coefficient at the finest scale 2"-'. A summarized 
description of our reconstruction procedure is given in the Algorithm[T] 

In the following we will use this algorithm to compute CSR reconstructions as well as A-CSR reconstructions, 
i.e., in the formulation of Algorithms [T] K may be the full or reduced system matrix. 

Eventually, we would like to note that the reconstruction Algorithm [T] is completely free of any parameter The 
weight sequence w appearing in the -penalized Tikhonov functional Q or ( |A-CSR[ ), respectively, is set adaptively 
during the iteration. 

7.2. Visibility of curvelets &- Dimensionality reduction 

In our first experiment, we are going illustrate the visibility of curvelets under the limited angle Radon transform 
■Ro (cf. Theorem |4.1| ) for different values of O. To this end we consider the function 

which is a linear combination of curvelets i/r,, / € { 1 , 2, 3, 4), at a fixed scale 2"** and orientations 9i - 0° , 62 = 20°, 63 - 
60° and 64 = 90°, see Figure|5] We computed the limited angle Radon transform H^sif - 'Rtt>i//i + 'R^ip2+ ^oiAs + ^o'A4 
and its inverse using the angular range parameters O = 35° and O = 80°. The results of this experiment are shown in 
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Algorithm 1 Reconstruction algorithm 



Require: < s< s„ < s <2/ \\Kf; 

1: 7 <— largest scale parameter of the curvelet decomposition; 

2: maxlter «— maximum number of iterations; 

3: c ^ 0; 

4: iter <- 0; 

5: while (iter < maxlter) do 

6: cr<- 1.4826 -MADCo); 

7: for each ( j, l,k) e I do 

8: Njj - number of curvelet coefficients at scale 2"-' and orientation 9jj; 

9: T«.w^230--^)/4^^21og,A^,;,; 

10: end for 

11: c^Sr:.{c-sK''iKc-y'^)); 

12: iter <— iter + 1; 

13: end while 



Figure |6] The first column shows the limited angle Radon transforms "Rqif of / for different values of O, whereas the 
second column shows the inverse Radon transforms 'R^'R,j, f from the corresponding limited angle data. In the first 
row we see that only those curvelets are visible in the reconstruction which correspond to 6i - 0° and 62 - 20°, i.e., 

5^35-^35"/ = ^1 + <A2. 

In the second row, we see that another curvelet 1^3 (corresponding to 63 - 60°) becomes visible by enlarging the 
angular range, i.e., 

^80-^80°/ = lAl +'A2 + 'A3. 

To explain this effect we computed the set of invisible curvelet coefficients according to ([TSj. As a result, we see that 

jmvisible ^ (3^4j ^jjj jmv.sible ^ 

As a rule of thumb, we can conclude that curvelets having orientations within the available angular range [-O, <!)] are 
visible for the angular range Radon transform. However, curvelets which correspond to missing directions are not 
visible for the limited angle Radon transform 



\\ 



Figure 5: A Matlab generated image of a function which is given as a linear combination of curvelets at scale 2 with orientations 
6 {0°,20°,60°,90°1. 
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Limited angle Radon transform 



Inverse Radon transform 



i 



(D = 35° 



W 



O = 80° 



\\ 



Figure 6: The first column shows the limited angle Radon transform (sinogram) of the image given in Figure|5]at an angular range 
[-O, O] for O = 35° and O = 80°. The second column shows the corresponding inverse Radon transform, i.e., 'R^'R^f. This figure 
illustrates the Theorem |4. 1| We can observe that the curvelet with the orientation 6 = 90° is invisible in both cases. Whereas, the 
curvelet with orientation 6 = 60° is invisible for O = 35° (first row), i.e., lies in the kernel of the limited angle Radon transform, 
but visible for O = 80° (second row). 



Now it is obvious that if an arbitrary function / is represented in terms of curvelet coefficients, we can seperate 
the visible and invisible parts of this function by 

/ = ^ C„(/r„ + ^ c„lf/„. 

„gjviabk „gjin,isible 

This separation depends only on the parameter The adapted dimension of the limited angle problem is then given 
by number of visible curvelets In the next experiment we computed the full and the reduced dimension for 

an image / of size 256 x 256 using CurveLab version 2.1.2, |3|. The results of this experiment are plotted in Figure 
|7] The dimension of the non-adapted problem in the curvelet domain is constant for all angular ranges. However, 
the dimension of the adapted problem shows a strong dependence on the available angular range. We can observe a 
significant dimensionality reduction for any angular parameter satisfying O < 150°. 

Moreover, we can observe a piecewise constant behavior of the reduced dimension. The dimension increases 
stepwise linearly as the angular range increases. The reason for this stepwise structure lies in the fact that curvelets 
remain visible as long as suppi^y / ;^ <^ ^ 0> see also Figure|4] The length of one such step therefore corresponds to 
the length of the of the support of the angular window V of curvelets at the finest scale 2 ■', i.e., to jsupp V ^^^'^^ ' -jj. 
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full dim 

I reduced dim 

■ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ >• 

0° 20° 40° 60° 80° 100° 120° 140° 160° 180° 

Angular range 

Figure 7: Dimension of the full problem ([32j, | | and of the adapted problem ([33j, | [ for an image of size 256 x 256. 

The plot shows the dependence of the full and reduced dimension on the available angular range [0, 0]. Since the full problem is 
formulated in terms of all curvelet coefficients, its dimension is given by the number of all curvelet coefficients The adapted 
problem, however, is formulated only in terms of visible curvelet coefficients. Hence, the reduced dimension, given by 
depends strongly on the available angular range. 



7.3. CSR vs. A-CSR: Execution times &■ reconstruction quality 




(a) Shepp-Logan head phantom 



(b) Brainstem | 
Figure 8: Original images. 



(c) Radial pattern 



In the following experiments we are concerned with limited angle reconstructions obtained via the adapted and 
via the non-adapted curvelet sparse regularization. In particular, we are going to investigate these reconstructions in 
terms of execution time of the reconstruction procedure and the reconstruction quality. 

Experimental setup 

The limited angle Radon transform was computed for test images which are shown in Figure [8] To this end, 
we considered different angular ranges [0,0], where the parameter was chosen to vary between 1° and 180°, 
i.e., G {1°, . . . , 180°). The generation of the limited angle data was done using the Matlab function radon. To 
simulate practical conditions, the generated data was corrupted by a white Gaussian noise, which was generated by 
the Matlab function randn. Having generated the limited angle data = '^o/ + we computed the CSR and A- 
CSR reconstructions using 100 iterations of the Algorithm [T] Instead of computing the system matrix directly and 
storing it in the memory, we implemented the transform Kc - 'R<s,T*c, ( [32] i, and its adapted version ( [33] l using the 
Matlab function radon and the CurveLab version 2.1.2, |3 |. Furthermore, we computed filtered backprojection (FBP) 
reconstructions using the Matlab function iradon. 
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Execution times 

We start by comparing the execution times of the CSR reconstructions to those of A-CSR reconstructions. The 
resuhs of this experiment are plotted in Figu re [9] h i this plot, the dotted line[ 
CSR reconstructions, whereas the solid line 



jindicates the execution times of the 

shows the execution times of the adapted approach ( |A-CSR| l. The 
dependence of the execution times on the available angular range shows in both cases a linear behavior. In particular, 
we can observe a significant speedup of the adapted procedure, especially for angular ranges [0, 0] with < 120°. 
The speedup exhibits a linear dependence on the available angular range which is due to the dimensionality reduction 
which was presented in the previous experiment, cf. Figure [7] 



600 



•3 400 



g 
■a 

3 



200 




90° 120° 
Angular range 

Figure 9: Execution times for CSR and A-CSR reconstruction using 100 iterations of Algorithm[T| Reconstruction of the Shepp- 
Logan head phantom of size 256 x 256 (Figure [^[a)[ l at different angular ranges [0, 0], 6 {1°, . . . , 180°). There is a significant 
speedup of the reconstruction procedure when using A-CSR. 



Reconstruction quality 



The results of the limited angle reconstruction are show in Figures 10-12 for angular ranges [0, 35°] and [0, 160°]. 



The original images corresponding to these reconstruction are shown in Figure |8] We investigate the reconstruction 
quality, first, by considering the CSR and the A-CSR reconstructions of the Shepp-Logan head phantom. These 



reconstructions are shown in the first and in the second column of Figure 10 By visually inspecting the images in 
each row separately, we can observe no difference in image quality. Inspecting reconstructions of the brainstem image 



(Figure 8b i and the radial pattern image (Figure 8c i which are shown in Figures 1 1 and 12 we can again observe that 
there no difference in image quality between CSR and A-CSR reconstructions. Therefore, we infer that the CSR and 
the A-CSR produce reconstructions of the same visual quality. 

To make these observations independent of visual perception, we used the mean squared error (MSE) as a quality 
measure. This is defined as 



MSE(c' 



1 ™ 
' N ^ 



where c are the curvelet coefficients of the original image and c' 

obtained via CSR or A-CSR at different angular ranges. The resulting MSE values are plotted in Figure 13 



"^^^ denotes those curvelet coeflicients which were 

As 



a function of the angular range parameter <b, MSE is decreasing for the non-adapted as well as for the adapted 
reconstruction method. However, the plots of the MSE values for CSR and A-CSR reconstructions again seem to be 



identical, cf. Figure 13 To refine our investigation we additionally consider the relative MSE, 



MSE(cCSR^^A-CSR) 



j_ Y I CSR 



A-CSR I- 



which compares the reconstructed curvelet coefficients obtained via CSR and those obtained via A-CSR. The plot 



of these values is shown in Figure 14 Here, we can observe how large the difference between the CSR and A-CSR 
reconstructions is in the case of Shepp-Logan head phantom reconstructions. Depending on the available angular 
range, the relative MSE values differ between 10"^ and 10"^. 
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As a result of the above discussion, we can conclude that the difference in the reconstruction quality of the CSR 
and the A-CSR reconstructions is very small. Visually, the reconstructions are not distinguishable. Therefore, the 
advantage of using the A-CSR approach consists in its significantly faster execution time. 

However, one might still ask where these differences come from? A possible explanation would be as follows: 
Since the reconstructed sequence of curvelet coefficient c'^^^ contains invisible curvelet coefficients, these values may 
be not zero after a finite number of iterations and, hence, these values would contribute to the relative MSE. Such a 
behavior was observed during these numerical experiments, thought the values of the invisible curvelets were very 
small. 



CSR 



A-CSR 



FBP 



= 35° 



= 160° 








Figure 10: Reconstruction of the Shepp-Logan head phantom of size 256 x 256 (Figure i a) i at an angular range [0, 0] and 
noiselevel 2% by using CSR, A-CSR and FBR In the above matrix of images, each row shows a reconstruction corresponding to 
the angular range parameter e {35°, 160°). Visually, there is no difference between CSR and A-CSR reconstructions for any of 
the angular ranges. However, the A-CSR reconstructions were computed significantly faster In contrast to FBP, CSR and A-CSR 
reconstructions appear less noisy. Though CSR and A-CSR reconstructions are slightly smoother than the FBP reconstructions, 
the edges are clearly visible. 



21 



CSR A - CSR FBP 


= 35° 


i 







= 160° 




Figure 11: Reconstruction of a brainstem glioma of size 301 x 310 (Figure i b) i at an angular range [0, 0], noiselevel 2% using 
CSR, A-CSR and FBR Again, there is no visible difference between CSR and A-CSR reconstructions for any of the angular ranges. 
Though CSR and A-CSR reconstructions are slightly smoother than the FBP reconstructions, the overall image quality is better in 
the case of CSR and A-CSR. 
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Figure 12: Reconstruction of a radial pattern of size 256 x 256 (Figure |q[c) [( at an angular range [0, 0], noiselevel 2% using CSR, 
A-CSR and FBR Same conclusions can be drawn as in the Figures [TO]and|lT] However, in this case we can additionally observe 
that only those lines are reconstructed whose normal directions are within the available angular range. 
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Figure 13: Mean squared error (MSB) of reconstructed curvelet coefficients, i.e., jj 2!l=i ~ '^'IT'^I ' where c denotes the vector 
of curvelet coefficients of the original image and c^"^ denotes the vector of reconstructed curvelet coefficients. This plot shows the 
results of the reconstruction of Shepp-Logan head phantom at an angular range [0, 0] and noiselevel 2%. 
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Figure 14: Relative MSE of reconstructed curvelet coefficients, i.e., jj 2"=i - c^^'^^'^p, where c'^^^ and c-^-csr ^jgjjQjg (j^g 
curvelet coefficient vector of the CSR reconstruction and A-CSR reconstruction, respectively. This plot shows the results of the 
reconstruction of Shepp-Logan head phantom at an angular range [0, 0] and noiselevel 2%. 



Eventually, we compare the reconstruction quality of CSR and A-CSR reconstruction to the quality of filtered 
backprojection (FBP) reconstructions. From Figures [TO] - [12] we can observe that the FBP reconstructions contain 
more noise than reconstructions obtained through curvelet sparse regularization. Visually, the FBP reconstructions 
seem to be inferior to the CSR and A-CSR reconstructions. On the other hand, the visual impression of the CSR and 
A-CSR reconstruction appears to be quite good. Though CSR and A-CSR reconstructions are slightly smoother than 
the FBP reconstructions, all details are well preserved and the edges are still clearly visible. 

To veiify the visual impressions, we computed the peak signal-to-noise-ratio (PSNR) of the normalized recon- 
struction^by 

PSNR(c'^'='=) = lOlogl — — 

These values are shown in the Table[T] For each test image and each angular range, we can observe that the PSNR val- 
ues of curvelet sparse regularizations (CSR and A-CSR) are considerably larger than those of the FBP reconstructions. 
Since larger PSNR values correspond to a better image quality, these results again confirm the visual impression. 
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(a) Shepp-Logan head phantom (b) Brainstem 
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Table 1: PSNR values of normalized reconstructions. 

7.4. Comments 

Our intention to perform these experiments was to give a practical proof of concept for our results. The implemen- 
tation of the reconstruction algorithms is therefore very rudimental and, hence, there is much room for improvements 



"The gray values of the reconstructed images were normalized to the interval [0, 1] 
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or optimizations. For example, the execution times that are presented in Figure|9]may be improved by a more elaborate 
implementation of the Algorithm[T] Though, there are many other algorithms available in the literature, the reason to 
use the iterative soft-thresholding algorithm for our experiments was its simplicity. 

8. Summary & Concluding remarks 

In this work we have introduced curvelet sparse regularization as a stable reconstruction method for the limited 
angle tomography. The stabilizing nature of the this method was demonstrated in numerical experiments. In com- 
parison to the FBP reconstructions, curvelet sparse regularization reconstructions offered a superior reconstruction 
quality. Another issue, that was addressed by CSR is its ability to produce edge-preserving reconstructions. The 
numerical experiments confirmed that to some extent. Wa have seen that all details in the in the CSR reconstructions 
were well preserved and the edges were clearly visible. However, CSR reconstructions were found to be smoother 
(more blurry) than FBP reconstructions. We believe that an even better edge-preservation can be achieved by tuning 
the reconstruction procedure. 

The main part of this work was devoted to the characterization of curvelet sparse regularizations in limited angle 
tomography. In Section |4j we have given a characterization of limited angle CSR reconstructions in terms of visible 
and invisible curvelet coefficients. Based on this characterization, an adapted CSR method was formulated. The adap- 
tivity of this approach results from the fact that, depending on the available angular range, the curvelet dictionary can 
be partitioned into a sub-dictionary of visible curvelets and a sub-dictionary of invisible curvelets. So, by formulating 
the reconstruction problem only with respect to the visible curvelet sub-dictionary, the problem becomes adapted to 
the limited angle geometry. Moreover, this entails a significant dimensionality reduction of the original reconstruction 
problem. This dimensionality reduction can be easily implemented in practice. A proof of this concept was given 
by numerical experiments. As a result, we found that the achieved dimensionality reduction is considerable, espe- 
cially, when the available angular range is small. Consequently, a significant speedup of the reconstruction algorithms 
was observed. The reconstruction quality of the adapted approach, however, was found to be equal to that of the 
non-adapted method. 

Furthermore, we would like to note that the results of this work can be generalized to the three-dimensional setting. 
The ideas of this work carry over to this situation, even though, the analysis is more technical in this case. 

We conclude this article by emphasizing the role of curvelets in case of limited angle tomography and summarize 
the reasons why they were used in this work: On the one hand, curvelets provide a sparse representation of functions 
with an optimal encoding of edges. These properties qualify curvelets for the use in sparse regularization and give 
rise to an edge-preserving reconstruction. On the other hand, curvelets are highly directional. Therefore, they enable 
a separation of visible and invisible structures of a function which is imaged at a limited angular range. Because of 
this directionality, curvelets allow to adapt the problem the limited angle setting. 

Acknowledgements 

The author gratefully acknowledges the support from GE Healthcare, Image Diagnost International, Munich. He 
especially thanks Peter Heinlein (GE Healthcare, Image Diagnost International, Munich) for his support during this 
work. The author also acknowledges the support of the TUM Graduate Schools Thematic Graduate Center ISAM at 
Technische Universitat Miinchen, Germany. 

References 

[1] Bredies, K., Kunisch, K.. Pock, T., 2010. Total Generalized Variation. SIAM Journal on Imaging Sciences 3 (3), 492 - 526. 
[2] Bredies, K., Lorenz, D. A., 2008. Linear Convergence of Iterative Soft-Thresholding. Journal of Fourier Analysis and Applications 14 (5-6), 
813-837. 

[3] Candes, E., Demanet, L., Donoho, D. L., Ying, L., 2008. Curvelab-2.1.2. http : //www. curvelet . org/. 

[4] Candes, E. J., Donoho, D. L., 2002. Recovering edges in ill-posed inverse problems: optimality of curvelet frames. Ann. Statist. 30 (3), 

784—842, dedicated to the memory of Lucien Le Cam. 
[5] Candes, E. J., Donoho, D. L., 2004. New tight frames of curvelets and optimal representations of objects with piecewise C~ singularities. 

Comm. Pure Appl. Math. 57 (2), 219-266. 



25 



[6] Candes, E. J., Donoho, D. L., 2005. Continuous curvelet transform. I. Resolution of tlie wavefront set. Appl. Comput. Harmon. Anal. 19, 
162-197. 

[7] Candes, E. J., Donoho, D. L., 2005. Continuous curvelet transform. II. Discretization and Frames. Appl. Comput. Harmon. Anal. 19 (2), 
198-222. 

[8] Caselles, V., Chambolle, A., Novaga, M., 2007. The discontinuity set of solutions of the TV denoising problem and some extensions. 

Multiscale Modeling & Simulation 6 (3), 879-894. 
[9] Chui, C. K., 1992. An introduction to wavelets. Vol. 1 of Wavelet Analysis and its Applications. Academic Press Inc., Boston, MA. 
[10] Daubechies, I., et al., 2004. An iterative thresholding algorithm for linear inverse problems with a sparsity constraint. Comm. Pure Appl. 

Math. 57 (II), I4I3-I457. 

[1 1] Davison, M. E., 1983. The ill-conditioned nature of the limited angle tomography problem. SIAM Journal on Applied Mathematics 43 (2), 
428-448. 

[12] Engl, H. W., Hanke, M., Neubauer, A., 1996. Regularization of inverse problems. Vol. 375 of Mathematics and its Applications. Kluwer 
Academic Publishers Group, Dordrecht. 

[13] Fadili, J. M., Peyre, G., 201 1. Total Variation Projection With First Order Schemes. IEEE Transactions on Image Processing 20 (3), 657-669. 

[14] Frikel, J., April 2010. A new framework for sparse regularization in limited angle x-ray tomography. Biomedical Imaging: From Nano to 
Macro, 2010 IEEE International Symposium on, 824—827. 

[15] Frikel, J., May 201 1. Short communication: Dimensionality reduction of curvelet sparse regularizations in limited angle tomography. Sub- 
mitted to the Proceedings in Applied Mathematics and Mechanics (2 pages). 

[16] Griesse, R., Lorenz, D. A., 2008. A semismooth newton method for tikhonov functionals with sparsity constraints. Inverse Problems 24 (3), 

035007 (19pp). 

URL http : //stacks . iop . org70266-5611/24/035007 

[17] Hansen, P. C, Sidky, E. Y., Pan, X., may 201 1. Accelerated gradient methods for total-variation-based CT image reconstruction. arXiv.org 
math.NA. 

[18] Herman, G. T., Davidi, R., 2008. Image reconstruction from a small number of projections. Inverse Problems 24 (4), 04501 1. 

URL http : //stacks . iop . org/0266-5611/24/i=4/a=045011 
[19] Kolehmainen, V, et al., 2003. Statistical inversion for medical x-ray tomography with view radiographs: II. Application to dental radiology. 

Phys. Med. Biol. 48, I465-I490. 

[20] Lorenz, D. A., Trede, D., 2008. Optimal convergence rates for tikhonov regularization in besov scales. Inverse Problems 24 (5), 055010. 

URL http : //stacks . iop . org/0266-5611/24/i=5/a=055010 
[21] Mallat, S., 2009. A wavelet tour of signal processing, 3rd Edition. Elsevier/Academic Press, Amsterdam, the sparse way. With contributions 

from Gabriel Peyre. 

[22] Natterer, R, 1986. The mathematics of computerized tomography. B. G. Teubner, Stuttgart. 

[23] Quinto, E. T, 1993. Singularities of the X-ray transform and hmited data tomography in K- andR'. SIAM J. Math. Anal. 24(5), I2I5-I225. 

[24] Quinto, E. T, 2006. An introduction to X-ray tomography and Radon transforms. In: The Radon transform, inverse problems, and tomogra- 
phy. Vol. 63 of Proc. Sympos. Appl. Math. Amer. Math. Soc, Providence, RI, pp. 1-23. 

[25] Radiopedia.org, 2010. 

URL http : //radiopaedii Torg/cases/bralnstem-gliomal 

[26] Rantala, M., et al., February 2006. Wavelet-based reconstruction for limited angle x-ray tomography. IEEE Transactions on Medical Imaging 
25(2), 210-217. 

[27] Ring, W., 2000. Structural properties of solutions to total variation regularization problems. Mathematical modelling and numerical analysis 
34 (4), 799-810. 

[28] Rockafellar, R. T., 1970. Convex analysis. Princeton Mathematical Series, No. 28. Princeton University Press, Princeton, N.J. 
[29] Scherzer, O., Grasmair, M., Grossauer, H., Haltmeier, M., Lenzen, F, 2009. Variational methods in imaging. Vol. 167 of Applied Mathemat- 
ical Sciences. Springer, New York. 

[30] Stein, E. M., Weiss, G., 1971. Introduction to Fourier analysis on Euclidean spaces. Princeton University Press, Princeton, N.J., princeton 
Mathematical Series, No. 32. 



26 



